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I. SUMMARY 


The multidimensional, ensemble-averaged, compressible, time-dependent Navier- 
Stokes equations have been used to study the turbulent flow field in two- and 
three-dimensional turbine cascades. The viscous regions of the flow were re- 
solved and no-slip boundary conditions were used on solid surfaces. The calcu- 
lations were performed in a constructive 0-type grid, which allows representa- 
tion of the blade rounded trailing edge. Converged solutions were obtained in 
relatively few time steps (-80-150) and comparisons of both surface pressure 
and heat transfer showed good agreement with the data. The three-dimensional 
turbine cascade calculation showed many of the expected flow field features. 



II. INTRODUCTION 


An important factor in the overall design of advanced gas turbine engines is 
the flow in and about the turbine blade passages. These components are a 
source of aerodynamic loads and losses that can control the overall machine 
efficiency and that may be subject to extreme heat transfer problems. Inaccur- 
ate estimates of both loss coefficients and heat transfer rates may lead to 
erroneous predictions of engine performance and to poor component life or even 
catastrophic failure of component parts, respectively. Therefore, an accurate 
estimate of both the flow field and the accompanying heat transfer properties 
in the turbine passages would be a valuable tool for the gas turbine design 
engineer. 

In recent years, reliable and efficient computational procedures for predicting 
the flow field and the accompanying heat transfer characteristics within the 
turbine blade passages have been developed. Several methods for analyzing the 
two-dimensional (2-D) passages in the turbine blades are available. Many 
analyses are based on solving inviscid flow equations (Ref 1 and 2) and, in 
many cases, these analyses are capable of predicting the blade pressure distri- 
bution well. 

Due to a lack of any viscous phenomena, however, these analyses cannot predict 
heat transfer or loss information. To predict heat transfer and/or viscous 
loss data, a set of boundary layer equations can be solved in conjunction with 
the inviscid analysis. A review of boundary layer techniques for turbo- 
machinery applications has been presented by McDonald (Ref 3). When the in- 
viscid flow boundary layer analysis technique is used, the calculations can be 
made in an interactive or noninteractive mode. In the noninteractive mode, an 
inviscid flow is calculated and the boundary layer is then calculated subject 
to the inviscid flow pressure distribution. This approach is viable if viscous 
displacement effects are small. When viscous displacement effects are signifi- 
cant, then the calculation should be performed in an interactive mode (Ref 4), 
which recognizes the mutual dependence of the inviscid and viscous solutions. 

An alternative method of determining the 2-D viscous flow in the turbine pas- 
sages is to solve an ensemble-averaged set of Navier-Stokes (N-S) equations in 
conjunction with a turbulence model to predict the entire passage flow field. 
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The N-S procedure has distinct advantages over other methods. The N-S analysis 
precludes any need to divide the flow into inviscid and viscous subregions. 
Although in some cases such a division may not be troublesome, in other cases 
no clear cut division is obvious. The N-S analysis solves the entire flow 
field with a single set of equations and, therefore, the influence of the 
viscous phenomena, taken in the boundary layer and in the wake, on the pressure 
distribution is inherent in the solution. Similarly, streamwise flow separa- 
tion is not difficult to solve with the N-S analysis. This fact is particu- 
larly important in cases involving distributed surface injection or blades sub- 
jected to high loading where local regions of separation, which may be regions 
of high heat transfer, may occur. Although boundary layer analyses can proceed 
through separation if approximations are made (Ref 5), such analyses contain 
serious approximations in and near the separated region if these regions are 
not small and thin. The results in this instance must be viewed with caution. 
The N-S analysis contains no approximations, other than those involved in 
turbulence modeling in this region. 

When considering three-dimensional (3-D) cascades, which include endwall and 
corner effects, the flow phenomena become even more complicated. In addition 
to possible 2-D flow-type situations in the vicinity of the midplane, complex 
3-D, viscous flow patterns appear in the endwall regions. These flow patterns 
include the leading edge horseshoe vortex at the strut endwall intersection, 
the complex vortex patterns passing through the passage (Ref 6), and the corner 
regions. These complex 3-D, viscous flow regions argue strongly for a full 
N-S approach. 

In simple quasi-two-dimensional situations, the viscous and turbulent near wall 
region velocity profile can sometimes be described by 3-D extensions of the 
law of the wall for turbulent boundary layers, although these 3-D extensions 
are much less universal than their 2-D counterparts. In the complex 3-D flows 
existing in the turbine passage, the near wall flow is not well described by 
such a wall law. Consequently, it becomes necessary to define the near wall 
region with mesh points when heat transfer and loss mechanisms are of interest. 
The added mesh point density and the presence of these additional length scales 
in the solution place very stringent requirements on the numerical methodology. 
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Successful calculations that applied a N-S procedure to several 2-0 cascade 
problems using no-slip boundary conditions (as opposed to a wall law) have been 
reported in Ref 7 through 12. The procedure used in these efforts was the con- 
sistently split linearized block implicit (LBI) scheme of Briley and McDonald 
(Ref 13). The numerical scheme is embodied in a general computer code termed 
MINT (multidimensional, implicit, nonlinear, time-dependent). The particular 
form of the code being used for the cascade application solves the general 
tensor form of the N-S equations and, therefore, can be used with a general 
coordinate system. The dependent variables in the analysis are the Cartesian 
velocity components, the density, and, if turbulent flow with a turbulence 
energy model is considered, the turbulence kinetic energy, k, and if a two- 
equation model is used, the dissipation rate, c. The choice of Cartesian 
velocity components as dependent variables is based on experience obtained with 
a variety of choices including the contravariant velocity components, the 
physical velocity components in the coordinate directions, and the Cartesian 
velocity components as discussed in Ref 14. 

To date, the cascade analysis has been used to predict the 2-0 flow in several 
cascade configurations using a C-type coordinate system. In Ref 10, flow 
fields were obtained for cascades formed by unstaggered NACA 0012 airfoils in 
subsonic laminar and turbulent flow. The results were obtained very economi- 
cally in relatively few time steps (~60 for laminar flow) and showed the ex- 
pected cascade flow field features. In Ref 11, the analysis was extended to 
the transonic regime and results were obtained for transonic flow through a 
Jose Sanz diffusion cascade. Once again, the results were obtained very eco- 
nomically in a relatively few number of time steps (-150) and showed the 
qualitative features expected in a transonic cascade flow field. 

In particular, the analysis showed its capability to predict transonic flows 
with sharp shock definition and, on a grid point-to-point basis, to cost little 
more than an efficient solution of the Euler equations. One case study cor- 
responded to conditions for which an inviscid, i.e. Euler equation, solution 
existed. The results of the N-S analysis showed good agreement with the in- 
viscid result. The differences observed were consistent with the inclusion of 
viscous boundary layer development in the N-S analysis. In comparing the vis- 
cous and inviscid results, with this particular N-S solution precedure no Kutta 
condition or exit flow angle specification is required, the exit flow follows 
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the specifications of the inlet total pressure and the exit static pressure. 
Inviscid flow calculations with other boundary conditions have the possibility 
of being different by virtue of these boundary conditions. Ref 7 and 9 contain 
calculations for both subsonic and transonic compressor configurations and a 
transonic turbine configuration. The predicted blade pressure distributions 
were comparable with the measurements of Stephens and Hobbs (Ref 15) and Hobbs, 
Wagner, Dennenhof f er, and Dring (Ref 16). The subsonic cascade flow (Ref 16) 
had two boundary layer profile measurements at the 97% chord pressure and suc- 
tion surface stations. Although only a tentative assessment of the procedure 
can be made on the basis of these two profiles in the immediate vicinity of 
the blade trailing edge, the comparison showed the predicted profiles to be in 
agreement with the measured data even though further turbulence modeling work 
with particular emphasis on transition was indicated. 

This effort applies this same N-S procedure to 2-D and 3-0 transonic turbine 
cascade flows. Such an application requires several new features, including a 
new coordinate capability, extension to three dimensions, and the inclusion of 
an energy equation in the governing system. In general, the geometrical con- 
figuration of the turbine blades impacts both the grid construction procedure 
and the implementation of the numerical algorithm. Because the turbine blades 
of interest, e.g. Turner and C3X turbine cascades, are characterized by very 
blunt leading edges, rounded trailing edges, and high stacking angles, a robust 
grid construction procedure that can accommodate the severe body shape while 
resolving regions of large flow gradients is required. A constructive 0-type 
grid generation technique that meets these requirements has been developed and 
used in this effort. An energy equation was activated in the code, 2-0 calcu- 
lations employing the N-S procedure were performed for the Turner and C3X tur- 
bine cascades, and the predicted pressure coefficients and heat transfer rates 
were compared with experimental data wherever available. A calculation with a 
film-cooling option applied over the C3X turbine blade surface was also demon- 
strated. In addition, the corresponding 3-D rectilinear C3X turbine cascade 
in which blade-endwall effects are present is considered. Throughout the ef- 
fort, no-slip boundary conditions with the viscous near wall definition of the 
flow were used as appropriate with the total pressure inflow and exit static 
pressure at outflow used to obtain the mass flow through the channel and to 
determine the flow turning. These boundary conditions were found the most 
plausible yet the least restrictive in terms of the predicted flow. 
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III. ANALYSIS 


This analysis is based on a solution of the ensemble-averaged Navier-Stokes 
(N-S) equations using the linearized block implicit (LBI) method of Briley and 
McDonald (Ref 13). The equations are solved in a constructive coordinate sys- 
tem (Ref 7) with density and the Cartesian velocity components as dependent 
variables. The application of the LBI method to the two-dimensional (2-D) cas- 
cade flow field problem together with the C-type grid coordinate system has 
been discussed in some detail in Ref 7 through 12. This effort deals with the 
application of the LBI method to the 2-D or 3-0 transonic turbine cascade flow 
field problem using an 0-type grid coordinate system and includes the effects 
of heat transfer. 

COORDINATE SYSTEM 

An important component of the turbine cascade analysis is the cascade coordin- 
ate system. Any coordinate system used in the analysis should satisfy condi- 
tions of generality, smoothness, resolvability, and allow easy application of 
boundary conditions. A coordinate system must be sufficiently general to allow 
application to a wide class of problems of interest if it is to be practical. 
The metric data associated with the coordinate system must be sufficiently 
smooth so that the variation from grid point to grid point does not lead to a 
numerical solution dominated by metric coefficient truncation error. This re- 
quirement differs from the requirement of the existence of a specified number 
of transformation derivatives. The coordinate system must resolve flow regions 
where rapid flow field changes occur. Finally, coordinates should allow accur- 
ate implementation of boundary conditions; for the cascade this requirement is 
equivalent to the requirement that the metric coefficients be continuous across 
the periodic lines where periodic boundary conditions are to be applied. 

To date, several types of coordinate systems are available. These include 
solutions based on a conformal transformation, solutions based on the solution 
of a Poisson equation (Ref 17), and constructive systems. This effort uses a 
constructive system based on the approach of Eiseman (Ref 18). Shamroth, et 
al (Ref 7 through 12) have applied the constructive technique to generate C- 
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type grids for a variety of cascades. The C-type grid requires an approxima- 
tion of the actual trailing edge geometry by means of a cusped geometry (see 
Figure 1). Therefore, the computed results based on the C-type grid represent 
only solutions for approximated geometry instead of realistic ones. The effect 
of the trailing edge approximation will depend on the type of cascade being 
considered. In compressor-type cascades, the trailing edge radius is small 
and the viscous boundary layers are subjected to large adverse pressure gradi- 
ents. Under these conditions, the blade viscous displacement effects are ex- 
pected to dominate any trailing edge geometry approximations and, consequently, 
the cusped approximation should not significantly influence the results for 
compressor blade calculations. The situation is not as clear, however, for 
turbine blades. Turbine blades have a large trailing edge radius and their 
boundary layers are subjected to strong favorable pressure gradients. In tran- 
sonic flow the pressure distribution and possibly the shock location are sen- 
sitive to small changes in the blade shape or the effective blade shape due to 
viscous effects. In these cases it is not clear if the trailing edge modifica- 
tion inherent in the C-type grid significantly affected results. As the first 
step in this effort, an 0-type grid capability was developed. 



Actual blade geometry 

Approximate blade geometry 

used for calculation 

TE84-058O 

Figure 1 . C-type grid for Turner cascade geometry. 
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A sketch of the 0-type grid coordinate system is presented in Figure 2. In 
brief, the coordinate system consists of a set of two families of curves; the 
l equals constant curves such as lines FG or H I in Figure 2 and the n 
equals constant curves such as ABCDEA or A'FB'HA' in Figure 2. The coordinate 
system is constructed by first forming the inner loop A'FB'HA 1 , which includes 
the blade. The blade may be specified either by an analytic equation or by 
discrete data points. If an equation is used then construction of the inner 
loop is straightforward. If the blade is specified by discrete data points 
then, in general, the points required on the inner loop will not coincide with 
any point used for blade specification. In this case, a curve fit is used to 
obtain the required inner loop points. The curve fit is based on a local para- 
bolic fit. For any given point required on the inner loop, a parabola is 
fitted through three adjacent specifying points, two on the right and one on 
the left, with the axis of the parabolic normal to the chord line connecting 
the two outer points. A second parabola is then fitted through the two points 
on the left and the one on the right. The location of the required point is 
obtained by means of a weighted average of these curve fits. The weighting 
factor is determined by the distance from the required point to the center 
specifying point of each parabola. This calculation is followed by construct- 
ing an outer loop ABCOEA, which consists of periodic lines BC and OE, a frontal 
curve CO, and a rear cap EB. Both the inner and outer loops are then repre- 
sented by parametric curves x = x(s) and y = y(s) where the parameter varies 
from zero to unity. The present coordinate generation process uses a two-part 
transformation for the inner loop. First, x and y are expressed as a function 
of s', the physical distance along the curve. Then s' is normalized so that 
its range is between zero and unity. Second, a transformation based on Oh (Ref 
19) is applied to the inner loop surface so that higher resolution at the lead- 
ing and trailing edge of the blade are achieved. The Oh technique uses a 
truncated series of error functions and complimentary error functions as the 
basis of the transformation. The feature of the technique is that the loca- 
tions that require higher resolution and their correspoinding grid point num- 
bers are specified explicitly through input. After applying Oh's transforma- 
tion, the desired grid point locations are obtained simultaneously. 

The outer loop is then parameterized to relate points on the outer loop to cor- 
responding points on the inner loop. For two outer loop points that are 
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periodic, such as G and I in Figure 2, to maintain periodicity, it is necessary 
that S„ = 1 - S T . These values will ensure that coordinate points on the 
upper and lower periodic lines will be periodically aligned. 

The present parameteri zation applies Oh's technique on an iterative basis over 
surface ABCC' so that the grid point locations are obtained as desired. The 
grid point locations for the surface C'DEA are obtained by means of correspond- 
ence to those on the surface ABCC 1 again on an iterative basis using Oh's tech- 
nique. This technique eliminates the complexity of the multipart transforma- 
tion technique used previously (Ref 7 through 11), and the grid point locations 
are obtained directly through Oh's transformation. 

Following the contruction and determination of the grid point locations for 
the inner and outer loops, two intermediate loops are constructed as shown in 
Figure 3. The first intermediate loop is constructed around the blade surface 
with a normal distance, h^ , from the inner loop. A similar loop is con- 
structed inside the outer loop with a normal distance, from the outer 
loop. Points on these intermediate loops must then be associated with a param- 
eter s, 0 < S < 1. This association is accomplished as shown in Figure 3 by 
setting S A , = S A and Sg, = Sg. These four loops allow the construction 
of the pseudoradial lines such as GF of Figure 2 by means of the multi loop 
method originally developed by Eiseman (Ref 20). 


B 



Figure 2. Constructive 0-type coordinate system. 
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Figure 3. Four basic loops of the constructive coordinate system. 

The multiloop method requires introduction of a position vector, P(r,s), with 
components (x, y), which will represent the coordinate location of grid points. 
Based on the four-loop construction process, vectors P*(s) are defined with 
i = 1, 2, 3, 4. Each PT has a coordinate (x^, y^) associated with it at 
specific values of s through the previously described contruction process. A 
radial parameter, r . , is then introduced. This parameter is defined, see 
Figure 3, as the distance from the loop in question to the inner loop normal- 
ized by the distance from the outer loop to the inner loop, h 3> Thus, = 0, 

r^ = h-j/h^ r 3 = (h 3 - h 2 )/h 3> r^ = 1 . With the definition of these 
quantities, the general position vector, P(r, s), is related to the loop posi- 
tion vectors, P^s), P^(s), P^(s), and P^(s), by means of 

P(r,s) «= (l-r) 2 (l-a |r )^(s) + ( 0 ,+ 2)(l-r) 2 (s) 

+ r 2 [l-a 2 (l-r)]^(s) + (a 2 +2)r 2 (l-r)l^(s) (1) 


where 




° 2 " 3(l-r 2 )-| 


It should be noted that at r = 0, P(0,s) = P^s) and at r = 1 , P(l,s) = 
P 4 (s). Further, because at r = 0, 


( 2 ) 
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-fji- (0,s) - [PjtsJ-P.ts^Oj+Z) (3) 

and at r = 1 

Cl. s) ■ [p 4 (s) - P 3 (s)](a 2 + 2) (4) 

specification of the derivatives at the inner and outer boundaries is deter- 
mined by the parametric representation of intermediate loops 2 and 3. Thus, 
the four-loop method allows specification of the boundary point locations and 
coordinate angles at these boundaries. This method of construction assures 
that the grid is orthogonal at the inner and outer loop boundaries. 

After loops 2 and 3 are constructed to satisfy the coordinate angle at the 
boundary points, the grid is constructed as follows. If the grid is to contain 
M pseudoradial lines (such as line FG of Figure 2) and N pseudoazimuthal lines 
(such as line ABCDEA) , the values of the pseudoradial coordinate are r(i) = 
i /( N - 1), i = 0, 1, 2, . . . , N - 1 and the values of the pseudoazimuthal co- 
ordinates are s(j) = j/(M - 1), j = 0, 1, 2, . . ., M - 1 then the position 
vector, i.e., the grid locations (x, y) for each point in the grid, is given 
by Equation (1 ) . 

The preceding discussion has assumed a uniform spacing in the radial direction. 
If radial grid point concentration is desired, it is necessary to assume a 
radial distribution function, such as 


tonhD(l-r) 

R » 1 — — — 

tanhD 


(5) 
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which concentrates points in the wail region. Grid points are then chosen at 
r(i) = (i)/(N-l) and the analysis proceeds as outlined. 

Grids generated with this procedure are shown in Figures 4 and 5. Figure 4 
represents a turbine cascade corresponding to the cascade of Turner (Ref 21). 
The second cascade, shown in Figure 5, corresponds to the C3X turbine cascade 
(Ref 22). 

GOVERNING EQUATIONS 

The ensemble-averaged, time-dependent N-S equations used in this project can 
be written in vector form as continuity 

-jf + V7>U = 0 (6) 

momentum 

dpi) —— _ _ 

+ V • ( pUU ) = — VP + V- (u + tt t ) (7) 

at 



Figure 4. 0-type grid coordinate system for Turner cascade. 


12 




Figure 5. O-type grid coordinate system for C3X cascade. 

and energy 

+ V(/>Uh) = -V- (0 + 0 T ) + — + & + pc 
dt Df 


(8) 


where p is density, U is velocity, p is pressure, « is the molecular stress 
tensor, is the turbulent stress tensor, h is enthalpy, Q is the mean 
heat flux vector, "q T is the turbulent heat flux vector, $ is the mean flow 
dissipation rate and c is the turbulence energy dissipation rate. If the 
flow is assumed as a constant total temperature, the energy equation is re- 
placed by 


T = T + — — = constant (9) 

♦ 2C p 

where T t is the stagnation temperature, q is the magnitude of the velocity, 

and C is the specific heat at constant pressure. In a number of cases con- 
P 

sidered in this report, constant total temperature has been assumed constant. 
This assumption was made to reduce computer run time where the constant T t 
assumption was warranted. Cases in which this substitution has been made are 
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identified in the description of the results. A number of terms in Equations 
6 through 8 require definition. The stress tensor appearing in Equation 7 is 
defined as 


v 


2 

2/iID — ( 3 >-K b )V- 


Ul 


( 10 ) 


where is the bulk viscosity coefficient, I is the identity tensor, and 
ID is the deformation tensor, defined by 


D =7((VU) + (VU) T ) (11) 

In addition, the turbulent stress tensor has been modeled using an isotropic 
eddy viscosity so that 


-T _ _ 


p u' u' = 2^- t D - — (/^ T V u + pK ) I 


( 12 ) 


where k, the turbulent kinetic energy, and uj, the turbulent viscosity, 
are determined by a suitable turbulence model. 

Equation (8) contains a mean heat flux vector defined as 

0 = — #c V T (13) 

and a turbulent heat flux vector defined as 

0 T = — k t VT < 14 ) 

where * and are the mean and turbulent thermal conductivities, re- 
spectively. 

The mean flow dissipation term, i, which also appears in Equation (8), is 
defined as 


<J> = 2/xD:ID-(-| / x-K )(V U) J 


(15) 
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In the equation of state for a perfect gas 


P = pRT 

where R is the gas constant, the caloric equation of state 

e = C V T 

and the definition of static enthalpy 

h = C p T 


( 16 ) 


(17) 


(18) 


supplement the equations of motion. 

Finally the flow properties y, k, and are determined using the follow- 
ing constitutive relations. 


The molecular viscosity, y, is determined using Sutherland's law, 

h- (JL \ 3/2 t q +s i 

Mo = ^V T+S, 


(19) 


where S 1 = 100°K for air. 

The bulk viscosity is assumed to be zero 


K b =0 (20) 

and the thermal conductivity is determined by use of a relation similar to 
Sutherland's law namely 



(J_f* T 0 + S g 

V V T+s * 


where S^ = 194°K for air. 


( 21 ) 
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DEPENOENT VARIABLES AND COORDINATE TRANSFORMATION 


The governing equations. Equations (6) through (8), are written in general 
vector form. Prior to their application to specific problems, it is necessary 
to decide on both a set of dependent variables and a proper coordinate trans- 
formation. Based on previous investigations (Ref 8), the specific scalar mo- 
mentum equations to be solved are the x, y, and z Cartesian momentum equations 
The dependent variables chosen are the physical Cartesian velocities u, v, w, 
and the density p. If the energy equation is solved, enthalpy is added to 
the set of dependent variables. 

The governing equations are then transformed to a general coordinate system in 
which the general coordinates, y^ , are related to the Cartesian coordinates, 
x 1 , x 2> and x 3 , by 


y J = y J (x, » * 2 > v*); j = '» 2 » 3 (22 

T = t 

As implied by Equation 22, the general coordinate, y\ may be a function of 
both the Cartesian coordinates and time. This coordinate time dependence will 
have an implication in so far as the choice of governing equation form is con- 
cerned. 


The governing equations can be expressed in terms of the new independent vari- 
ables, y J , as 


aw 

dr 


„ 

: td£ 

dw 


dw 


. dF 
* d£ 

+ 

dG 

dZ 

c <*H 

dF 
^x drj 

+ 7 >y 

dG 

d l 

SH 

+ ’’l d v 

dF 


dG 

+ t 
^dC 

x 

H 


h[ 


r ap ' * 


d^ 

d l 


dF, 

K 




y d£ 


dG| dG| 

+ T )^ 


dr, 


dG| 




dH, 

Jt +7 ?z 


dH, 




dH, 

K 


( 23 ) 
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through a straightforward application of chain rule differentiation. In Equa- 
tion (23) 


t ■ y' 

Tj = y 2 

C = y 5 


and 


p 


pu 


■ py 

pu 

F = 

pu 2 + P 

G = 

p uv 

P v 


p uv 
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p w 
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_ J 


pvw 
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T xy 
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T yy 


T yz 
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The governing equations have been written in this form in Ref 7 through 12; 
the form is termed the semistrong conservation form in Ref 23, and the chain 
rule conservation form in Ref 24. The metric coefficients do not appear within 
the derivatives. The experiences of Ref 7 through 12 show that results ob- 
tained by means of the semistrong conservation form are less sensitive to the 
method used to evaluate the metric coefficients than are results obtained by 
means of strong conservation form. The semistrong conservation form was used 
in this effort. 
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TURBULENCE MODELING 


Several alternative turbulence models have been considered in the course of 
this effort. In general terms, the models used were zero-, one-, and two-equa- 
tion models. 


Zero-Equation Model — Mixing Length 


Of all of the available turbulence models, Prandtl's mixing length model is 
probably still the most widely used. The model was originally developed for 
use in unseparated boundary layer flow situations and has been shown to perform 
well under such conditions. An economical advantage of the method is that it 
does not require additional transport equations to model the effect of turbu- 
lence, but rather relates the Reynold's shear stress to mean flow quantities 
by means of 


7>V] =Mt 


auj 

ax; 


(24) 


where 

/z T = pl 2 (21D : D) ,/2 


where 


i = 


min 




where d is the normal distance to the nearest wall and D is the van Driest 
damping coefficient given by 


D 


1 - exp(- y + /A + ) 


l m - 0.096 
•c - 0.4 

y + * du T /v 

u T - (y P ) 1/2 


(25) 
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where the local shear stress t is obtained from 

T t = (2D: D) ,/2 (26) 

and ID is defined by Equation (11). 

One problem in the mixing length formulation is the definition of 6. In 
boundary layers the streamwise velocity, U, approaches an edge velocity, U e> 
asymptotically, however a monotonic approach to an asymptotic edge velocity is 
not characteristic of N-S solutions. To avoid the ambiguity of determining 
the boundary layer edge, 6, as defined in the usual boundary layer context, 
i.e., 6 is the distance from the wall at which U/U e = 0.99, the following 
relation is used 


8 = 


2.0 d 


(q/q =c) 
m max 


(27) 


In other words, 6, is taken as twice the distance (measured away from the 
nearest wall) for which q/q = c. The value of c used in the present ef- 
fort was 0.81 . 

The model used in the wake is also a mixing length model in which the mixing 
length was made proportional to the wake height, 6, and a linear growth 
of 6 with distance was assumed based on the classical free jet boundary re- 
sults (Ref 25). With the free jet boundary growth assumption 


8 = (S ps +S ss ) + (0 - 2)(x ' x te ) (28) 

where 6 p$ and 6 $s are the pressure and suction surface trailing edge 
boundary layer thickness and x is the trailing edge location. The mixing 
length. It, was taken as 0.26. 
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One-Equation Model— k-8. 


Although the mixing length concept is valid for a variety of flows, some im- 
portant flow situations arise in which a less restrictive model is required. 

One such model is the one-equation turbulence model (Ref 26) in which a trans- 
port equation for turbulence kinetic energy, k, is formulated 

^- + V-(dUK)s V-(^IVK) + 2/z t OD: [D)-/>€ (29) 

at c K 

where (following Ref 27) o k is set to 1.0, and k is the turbulence kinetic 
energy 

K ~~Z u'u’ 

and the Prandtl -Kolmogorov relation, Equation (30), defines the turbulent 
viscosity as 


. P K 2 
H sZ v-—r- 


(30) 


In addition, the turbulence dissipation rate e is determined by 


€ 


.3/4 K 3 ' 2 

* i 


(31) 


where l is a relevant turbulent length scale for the problem of interest. 

The k-i. model has an advantage over the mixing length scale model in that 
the use of a transport equation for turbulence kinetic energy allows for its 
convection, production, and dissipation. This is important because it allows 
a nonequilibrium effect on the turbulence to be included in the calculation 
while the mixing length model can only account for local equilibrium turbulence 
effects by means of its association with the mean velocity field. A major dis- 
advantage that the k-8. model shares with the mixing length model is the re- 
quirement of length scale specification. Typically the mixing length, as de- 
scribed previously, is used as the representative length scale. 
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There are two approaches to modeling the flow in the near wall region where 
low local turbulence Reynolds numbers occur. The first is the wall function 
approach that does not resolve the near wall region but assumes specific func- 
tion forms for the required turbulence quantities and uses these forms to 
create the required normal derivative formulations at the first grid point from 
the wall. Such an approach obviously requires a detailed knowledge of the 
turbulence model dependent variables in the vicinity of the wall. Although 
reasonable function formulations can be specified for simple 2-D flows such as 
constant pressure boundary layers, specification in the much more complex flows 
is more difficult. Therefore, the alternative approach in which the viscous 
sublayer is resolved has been used. The method makes no approximation at the 
boundary, but requires that the near wall low turbulence Reynolds number 
physics be modeled. In this effort, a near wall model, which was successfully 
used by Shamroth and Gibeling (Ref 28) in a time-dependent airfoil flow field 
analysis, has been implemented in the computer code. The analysis of Ref 28 
follows the integral turbulence energy procedure of Ref 26 by using a turbu- 
lence function, a^ , where 


a 


1 


1 1/2 
2 C p 


(32) 


and a 1 is taken as a function of a turbulence Reynolds number of the form 


■l ‘ *o [t 5S^]/[ 10 * 6 66 *o - *)] 

a - .0115 
o 

where 

o 22 

f(R x ) - 100 R t °‘ ZZ 
f(R T ) - 68.1 R t ♦ 614.3 

and a cubic curve was fitted for values of R^ between 1 and 40. In this 
effort, R was defined as the local ratio of turbulent to laminar viscosity, 
a^ was evaluated by means of Equation (33) and related to a^ by means of 
Equation (32). As previously discussed, Ref 26 used an integral form of the 
turbulent kinetic energy and, therefore, R^ was defined as an average value. 


R < 1 

T — 

R x >1 40 


(34) 
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Two-Equation Model — k -c 


Although the one-equation approach does relieve some restrictions in the mixing 
length approach, it still requires specification of a length scale. The k-e, 
two-equation turbulence model (Ref 29 and 30), in which both the turbulence 
kinetic energy and the turbulence dissipation rate are governed by transport 
equations, represents a more general model. In this approach, the k-equation 
is as given in Equation (29), but the algebraic relation for e given by Equa- 
tion (31) is replaced by the following transport equation. 

dp€ +V-(/)U€)= V (J^IV«) + C,(2/i t ID: ID)!. +2fifi T (V 2 U) 2 -C z pl_ (36) 

at CT C K K 

However, attempts to solve Equations (29) and (36) without modification present 
problems because an appropriate boundary condition for c at a solid boundary 
is difficult to prescribe such that Equation (36) is satisfied. Following the 
suggestion of Jones and Launder (Ref 27), the turbulence dissipation equation 
has been modified by the inclusion of the term 

- 2/z^ t (V 2 U) 2 


in the energy dissipation equation, Equation (36), and by the inclusion of the 
term 


-2pv(VK u2 ) 2 


in the turbulence energy equation which then becomes 


dp K — ii— 

+ V-(/> UK) = V-(-^VK) + 2/i t (ID:ID)-/>€-2/>i/(VK ,/2 ) 2 < 37 ) 
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These additional terms allow an c = o wall boundary condition to be applied 
and appear to correctly model the near wall region as discussed in Ref 27. 

Following Ref 27, the following empirical relations are used 

a € c 1.3 
cr k = 1.0 
C ( = 1.43 

= 0.09 exp[-2.5/(l + R t /50)] 

C £ = 1.92 [l.O - 0.3 exp (~R*)] 

and R t is defined as 


P ** 



The turbulent eddy viscosity is evaluated by Equation (30). 

BOUNDARY CONDITIONS 

Boundary conditions play an important role in determining accurate solutions 
and rapid numerical convergence when solving N-S equations. The boundary con- 
ditions used in these calculations follow the suggestion of Briley and McDonald 
(Ref 31), which specifies upstream total pressure and downstream static pres- 
sure conditions. For the 2-D cascade system shown in Figure 2, BC and ED are 
periodic boundaries and periodic conditions are set here. 

Specification of upstream and downstream conditions is somewhat more difficult. 
For an isolated cascade, boundary conditions for the differential equations 
may be known at both upstream infinity and downstream infinity. However, since 
computation economics argues for placing grid points in the vicinity of the 
cascade and minimizing the number of grid points far from the cascade, the up- 
stream and downstream computational boundaries should be set as close to the 
cascade as is practical. In addition, with the particular body-fitted coordin- 
ates used, as the upstream boundary moves further upstream, the angle between 
pseudoradial and pseudoazimuthal coordinate lines becomes smaller. Decreasing 
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the coordinate angle causes the coordinate system to become less well-condi- 
tioned, increases truncation error (Ref 32), and increases the role of cross- 
derivative terms in the equations. All of these characteristics could be 
detrimental to the present numerical procedure and, therefore, they also argue 
for placing the upstream boundary as close to the cascade as possible. How- 
ever, when the upstream boundary is placed close to the cascade, most flow 
function conditions on the boundary will not be known, because these will have 
been changed from values at infinity by the presence of the cascade. 

In this effort, total pressure is set on boundary CC'D, shown in Figure 2 (Ref 
31). Unless boundary CC'D is very far upstream, the flow velocity along this 
boundary will not be equal to the velocity at upstream infinity because some 
inviscid deceleration will have occurred. However, as long as the boundary is 
upstream of the region of any significant viscous or shock phenomena, the total 
pressure on this boundary will be equal to the total pressure at upstream in- 
finity. Hence, total pressure is an appropriate boundary condition realisti- 
cally modeling the desired flow condition. In addition to specifying upstream 
total pressure, it is necessary to specify the inlet flow angle. In these cal- 
culations, a value was assumed constant on the upstream boundary at a specified 
value, although it is possible that in future studies the inflow angle distri- 
bution could be obtained from an inviscid calculation. The third condition 
set on the upstream boundary concerns the density. A zero density derivative 
at this boundary was specified as a numerical treatment of the boundary. 

The downstream boundary EAB is considered to be far away from the blade surface. 
A small portion of the boundary may contain inlet flow depending on the turning 
of the passage flow. For this inlet flow portion, the flow variables were set 
equal to those at the end of the periodic line, e.g. at point E in Figure 2. 

For outflow portion, the boundary was treated by setting a constant static 
pressure as a boundary condition, and by extrapolating (first derivatives) both 
velocity components along exit flow direction at this location. In the present 
application, a constant static pressure was set at downstream infinity, and 
hence it is assumed that the downstream boundary is located in a region where 
pressure is uniform. 
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Both the upstream and downstream boundaries have boundary conditions that are 
nonlinear functions of the dependent variables associated with them. These 
conditions are the specifications of total pressure on the upstream boundary 
and static pressure on the downstream boundary. These nonlinear boundary con- 
ditions are linearized in the same manner as the governing equations (by means 
of a Taylor expansion of the dependent variables in time), and then solved im- 
plicitly along with the interior point equations. No-slip conditions (except 
the film cooling option) together with zero pressure gradient were set at solid 
walls. If film cooling is applied on the blade surface, a wall velocity is 
specified by means of input for the portion of that surface. 

For a three-dimensional (3-D) rectilinear cascade configuration, which consists 
of a turbine cascade situated in the azimuthal -radial plane and bounded in the 
transverse direction by an endwall and a symmetry plane (see Figure 6), the 
incoming flow at the upstream boundary is in boundary layer form due to the 
endwall and not a uniform stream as in the 2-D unbounded cascade. A two-layer 
velocity profile condition in place of a uniform velocity profile condition 
was set at the upstream boundary. The so-called two-layer model (Ref 31) used 
at the inflow boundary is essentially a total pressure boundary condition ap- 
plied to the core flow with a specified boundary layer profile shape for the 
wall region. Matching the static pressure at the edge, defined by the first 
computational point from the wall at which |u|/|ul max was greater than or 
equal to 0.99 on the previous time step, enables calculation of I u I at this 
point. This calculation provides the required normalizing value for the pre- 
specified boundary layer profile shape. Overall, the method provides a mech- 
anism for drawing mass flow to satisfy the downstream pressure given an up- 
stream core total pressure while maintaining a given inlet endwall boundary 
layer shape. This specification corresponds to a wind tunnel experiment where 
stagnation conditions are set in an upstream reservoir and static pressure is 
set at some downstream location. No-slip conditions in conjunction with zero 
pressure gradient were set on the endwall. Symmetry boundary conditions were 
applied on the symmetry plane. The other necessary boundary conditions were 
treated the same as those for the 2-D case. 

The present formulation contains several important advantages over alternative 
formulations. Specification of upstream total pressure and downstream static 
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Figure 6. Three-dimensional C3X cascade with endwall 


pressure allows the flow to develop in a natural manner with no need to specify 
the velocity magnitude on the upstream boundary. Proper specification of the 
velocity magnitude on the upstream boundary may be difficult for 2-0 cases and 
is difficult in transonic cascades and/or 3-D cascades. In addition, specifi- 
cation of no-slip conditions at solid surfaces eliminates the problem of speci- 
fying a wall function. Although wall functions may be appropriate for rela- 
tively simple 2-D flow, their use becomes questionable for 2-D separated flows. 
In 3-0 flows, near wall flows are not well described by such universal laws 
and, consequently, their use in these cases does not seem appropriate, particu- 
larly if heat transfer or loss information is required. Use of the no-slip 
condition along with proper resolution of the boundary layer does not hinder 
convergence properties of the numerical method. 


NUMERICAL PROCEOURE 


The numerical procedure used to solve the governing equations is a consistently 
split LBI scheme originally developed by Briley and McDonald (Ref 13). A scheme 
similar in concept has been developed for 2-0 magnetohydrodynamic problems by 
Lindemuth and Killeen (Ref 33) as discussed in Ref 13 and 34. The method can 
be outlined briefly as follows: the governing equations are replaced by an 
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implicit time difference approximation, either a backward difference or Crank- 
Nicolson scheme. Terms involving nonlinearities at the implicit time level 
are linearized by Taylor expansion in time about the solution at the known time 
level, and spatial difference approximations are introduced. The result is a 
system of multidimensional coupled (but linear) difference equations for the 
dependent variables at the unknown or implicit time level. To solve these dif- 
ference equations, the Oouglas-Gunn (Ref 35) procedure for generating altern- 
ating direction implicit (ADI) schemes as perturbations of fundamental implicit 
difference schemes is introduced in its natural extension to systems of partial 
differential equations. This technique leads to systems of coupled linear dif- 
ference equations with narrow block-banded matrix structures, which can be 
solved efficiently by standard block-elimination methods. 

The method centers around the use of a formal linearization technique adapted 
for the integration of initial-value problems. The linearization technique, 
which requires an implicit solution procedure, permits the solution of coupled 
nonlinear equations in one space dimension (to the requisite degree of accur- 
acy) by a one-step noniterative scheme. Because no iteration is required to 
compute the solution for a single time step, and because only moderate effort 
is required for the solution of the implicit difference equations, the method 
is computationally efficient. This efficiency is retained for multidimensional 
problems by using block ADI techniques. The method is also economical in terms 
of computer storage, requiring only two time levels of storage for each depend- 
ent variable. In addition, the block ADI technique reduces multidimensional 
problems to sequences of calculations that are one-dimensional (1-D) in the 
sense that easily-solved narrow block-banded matrices associated with 1-0 rows 
of grid points are produced. A more detailed discussion of the solution pro- 
cedure as discussed by Briley, Buggeln, and McDonald (Ref 36) is given in Ap- 
pendix A. 

ARTIFICIAL DISSIPATION 

Due to frequent high Reynolds numbers typical of normal turbomachinery applica- 
tions, it is necessary to suppress spatial oscillations associated with central 
spatial difference approximations. Such suppression can be achieved by means 
of a dissipative spatial difference formulation (e.g., one-sided difference 
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approximations for first derivatives) or by explicitly adding an additional 
dissipative-type term. For the N-S equations, the latter approach has been 
selected because when an additional term is explicitly added, the physical ap- 
proximation being made is clearer than when dissipative mechanisms are con- 
tained within numerical truncation errors. Explicit addition of an artificial 
dissipation term also allows greater control over the amount of nonphysical 
dissipation being added. The most desirable technique would add only enough 
dissipative mechanism to suppress oscillations without deteriorating solution 
accuracy. Various methods of adding artificial dissipation were investigated 
in Ref 11, and these methods were evaluated in the context of a model 1-0 
problem containing a shock with a known analytic solution (1-D flow with heat 
transfer). The methods that were considered included second-order dissipation, 
fourth-order dissipation, and pressure dissipation techniques. 


As a result of this investigation, it was concluded that a second-order aniso- 
tropic artificial dissipation formulation suppressed spatial oscillations with- 
out adversely impacting accuracy and could be used to capture successfully the 
nearly normal shocks that are expected in transonic cascades. In this formula- 
tion the terms 

n-1 3 / 3<A n-1 3 / 3<t> 

p — (d — ]♦ p — jd — 

3x \ x 3x / 3y y y 3y 

are added to the governing equations where <t> = u, v, h, and p for the x- 
momentum, y-momentum, energy, and continuity equations, respectively. The ex- 
ponent n is zero for the continuity equation and unity for the momentum and 
energy equations. The dissipation coefficient d is determined as follows. 

The general equation has an x-direction convective term of the form a 3<t>/3x 
and an x-direction diffusion term of the form d(bd<t>/3x)/3x . The diffusive 
term is expanded 


d(b dcp/dxVdx 


bd^tfi/dx^ + db/dx d<fi/dx 


(38) 


and then a local cell 


Reynolds number, Re A x. is defined for the x-direction by 
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o - db/ dx Ax/b 


(39) 


Re 


Ax 


where b is the total or effective viscosity including both laminar and turbu- 
lent contributions and Ax is the grid spacing. The dissipation coefficient 
d x is not negative and is chosen as the larger of zero and the local quantity 
b (o Re - 1). The dissipation parameter is a specified constant and 
represents the inverse of the cell Reynolds number below which no artificial 
dissipation is added. The dissipation coefficient d y is evaluated in an 
analogous manner and is based on the local cell Reynolds number Re^y an d grid 
spacing Ay for the y-direction and the specified parameter o y . 

The question arises as to which values of a and a should be chosen. This 

x y 

question was assessed both through the model problem (Ref 11), and through cal- 
culations for a Jose Sanz compressor cascade (Ref 7 and 11). These results 
indicated that values of o = 0.5, which corresponds to a cell Reynolds number 
2 limitation, would severely damp physical variations. However, when a was 
set in the range 0.025 < a < 0.05, which corresponds to a cell Reynolds num- 
ber range between 40 and 20, spurious spatial oscillations were damped with no 
significant change in the calculated results as a was varied in this range. 

As discussed in Ref 7 through 11, the results showed good agreement with data. 
This agreement has since been confirmed at several other studies at Scientific 
Research Associates (SRA) such as 2-D and 3-D transonic nozzle flows (Ref 37) 
where a maximum acceptable value of a - 0.10 has been noted for most prob- 
lems. In cases where spatial resolution may be marginal, such as at the lead- 
ing edge of a relatively sharp edged blade, it may be necessary to increase 
a in this local area. However, a can be decreased to 0.10 or below if com- 
putational grid points are added in this region. 
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IV. TEST CASES AND RESULTS 


Several test cases were run to demonstrate and assess the computer code for 
the Navier-Stokes (N-S) analysis. The first case examined the two-dimensional 
(2-D) Turner turbine cascade (Ref 21). The second case examined the 2-D low 
solidity C3X turbine cascade (Ref 22). The final case considered was a three- 
dimensional (3-D) rectilinear C3X cascade. 

CASE 1— 2-D TURNER TURBINE CASCADE 

To assess the impact of the new 0-type grid on the calculated turbine cascade 
flow fields, a series of N-S computations were run for the Turner cascade with 
a rounded trailing edge. The Turner 0-type grid is shown in Figure 4 and con- 
tains no approximations in the trailing edge geometry. The results of these 
calculations were compared with the results obtained with the C-type grid for 
the Turner case with a cusped trailing edge (Figure 7) under the same flow con- 
ditions. As shown in Figure 1, the cusp is added by increasing the chord to 
keep as much of the geometry exact as practical. 

The C-type coordinate system, shown in Figure 7, consists of 30 points in the 

pseudoradial direction and 113 points in the pseudoazimuthal direction. The 

periodic boundary is extended 0.75 chords upstream of the leading edge. The 

0-type coordinate system, shown in Figure 4, consists of 30 points in the 

pseudoradial direction and 100 points in the pseudoazimuthal direction. The 

periodic boundary is extended 0.9 chords upstream of the leading edge and 1.0 

chords downstream of the trailing edge, respectively. The first coordinate 

-5 

line off the airfoil is placed at approximately 4 x 10 chords from the 
surface for both C- and 0-type coordinate systems, thus obtaining adequate 
boundary layer resolution. The geometric inlet angle for the Turner cascade 
is approximately 0 deg and the geometric exit angle is approximately 62 deg. 

The blade pitch to chord ratio is approximately 0.65 and the blade shape has 
a relativaely large radius of curvature at the leading edge. 

The first flow considered was assumed to be subsonic, laminar, and at constant 
total temperature. The calculation was run with a free-stream Reynolds number 
based on chord and inlet conditions of approximately 400 and a ratio of inlet 
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Figure 7. C-type coordinate system for Turner cascade. 
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total pressure to exit static pressure of approximately 1.25. The calculations 
were initiated by assuming uniform flow with a boundry layer correction on the 
airfoil surface. With these initial conditions, it took about 80 time steps 
to reach converged solutions for both C- and 0-type coordinate systems. In 
all cases, convergence was monitored by noting the decrease in maximum flow 
field residual. Based on previous experience, the solution can be considered 
converged when the maximum residual has decreased by three orders of magnitude 
from its initial value and only small changes are occurring in the calculated 
pressure and velocity distributions for a range of time step sizes. In Figure 
8, the pressure coefficient of the two coordinate systems is compared. The 
C-type grid blade is longer due to the cusped trailing edge. The results are 
in excellent agreement for both the suction and pressure surfaces except in 
the vicinity of the trailing edge, which is not unexpected in light of the dif- 
ferent trailing edge geometries for the two calculations. 



Figure 8. Pressure coefficients distribution of the subsonic laminar cascade. 
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The second flow considered was similar to the laminar one except the free- 
stream Reynolds number based on chord and inlet conditions was approximately 
0.527 x 10 6 and the flow was assumed to be turbulent. Using the convergent 
solution obtained from the laminar flow calculation as the initial condition 
and the mixing length turbulence model, 60 time steps were required to reach a 
converged solution. Good agreement was obtained for the turbulent subsonic 
case as can be seen in Figure 9. 

The final flow considered was a transonic turbulent flow corresponding to 
Turner's experiment (Ref 21). The free-stream Reynolds number based on chord 
and inlet conditions was approximately 0.884 x 10^ and the ratio of the inlet 
total pressure to the exit static pressure was approximately 1.778. Total 
temperature was assumed to be constant in the flow field. The calculations 
used the solution of the turbulent subsonic case as the initial condition, then 
the inlet total to exit static pressure ratio was increased linearly from 1.25 
to 1.778 within 15 time steps. Using the mixing length turbulence model, con- 
vergent solutions were obtained within about 80 time steps for both coordinate 



Figure 9. Pressure coefficients distribution of the subsonic 

turbulent cascade. 
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systems. In Figure 10, the surface pressure distribution (in terms of velocity 
ratio) for the two calculations is compared with the data of Turner (Ref 21). 
Both calculations are comparable with one another over the entire airfoil and 
both agree with the data over most of the surface except near the trailing 
edge. In view of the excellent agreement between the two calculations, it is 
concluded that the discrepancy between the measured data and C-type grid cal- 
culations originally noted in Ref 7 are not due to the approximated trailing 
edge geometry in the C-type grid calculation. In addition, the calculation 
confirms the operation of the N-S code in the 0-type grid mode. 

An inspection of Figure 10 reveals that both the 0- and C-type grid calculations 
underpredict the pressure in a similar manner in the region of 60-80% chord on 
the suction surface. There are several possible sources for this discrepancy. 
Computationally, the turbulence transition model could give rise to the ob- 
served differences. Experimentally, one possible source could be endwall ef- 
fects changing the axial velocity density ratio (AVDR) from unity. If this 
change were to occur, it would be particularly important at transonic speeds 
when the flow field becomes sensitive to the effective area ratio and neglect 
of a nonunity AV0R could lead to significant discrepancies. In view of the 
absence of endwall effect discussion in Ref 21, which focused primarily on 
transitional boundary layers, it is difficut to assess the actual source of 
the discrepancy. Because the intent of this effort was to verify the suit- 
ability of the 0-type grid, no further investigation of this question was un- 
dertaken. 

CASE 2 — 2-D C3X TURBINE CASCADE 

The coordinate system for the C3X turbine cascade is shown in Figure 5. This 
0-type grid consists of 30 points in the pseudoradial direction and 120 points 
in the pseudoazimuthal direction. In keeping with the objectives of the con- 
struction procedure, the minimum coordinate intersection angle is 20 deg and 
the upstream boundary is placed at 2.25 axial chords upstream of the leading 
edge and the downstream boundary is placed at 2.65 axial chords downstream of 
the trailing edge. High radial resolution is obtained near the surface of the 
blade, and the first coordinate line is located at a distance of 1.0 x 10 6 
axial chords from the surface. In addition, high pseudoazimuthal resolution 
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Figure 10. Turner cascade pressure distribution of the transonic case 


is obtained in both leading and trailing edge regions. The C3X cascade geome- 
try is given in detail in Ref 22. The geometric inlet angle is approximately 
0 deg and the geometric exit angle is approximately 72 deg. The vane spacing 
to axial chord ratio is approximately 1.5. The true chord to axial chord ratio 
is approximately 1.85. 

The first flow examined corresponded to case 143 in Ref 22 in which the inlet 
Mach number, Ma ] , was 0.17, the exit Mach number, Ma exit , was 1.05, the 
inlet Reynolds number based on true chord, Re 1 , was 0.63 x 10 & , and the 
estimated ratio of inlet total to exit static pressure, p t / p exit was 2 -°- 
The exit data given in Ref 22, i.e., the average Mach number downstream of the 
trailing edge, is not suitable for specifying boundary conditions at the down- 
stream boundary of the computational domain as required by the N-S procedure 
(2.65 axial chords downstream of the trailing edge). In this computational 
domain, the flow is subsonic at the rear cap. Thus a static exit pressure can 
be specified. However, a sensitivity study of the effect of p t / p ex ^ t on 
surface pressure distribution was undertaken, where P^ is the upstream stag- 
nation pressure, and P .. is the static pressure at the downstream boundary 

exit 

of the computational domain. Two values of p -j/ p ex -j£ were chosen, 1.9 and 
2.0. Constant total temperature was assumed in the calculation. The computed 
pressure distribution is shown in Figure 11 and is compared with the Allison 
experimental data (Ref 22) and the inviscid predictions (Ref 1) due to Oelaney. 

Figure 11 indicates that the pressure distributions on the pressure side of 

the turbine blade and the forward portion of the suction side are relatively 

insensitive to P./P ... However, the pressure distribution on the aft 
t exit 

portion of the suction side is sensitive to the pressure ratio, in that a 5% 
change in PVP ... results in a commensurate pressure variation on the 
blade's surface. Such behavior is not unexpected because the flow through the 
cascade is in the transonic regime and indicates a need for definitive specifi- 
cation of boundary conditions if a data comparison is to be made. 

Mixing length and two-equation k -c turbulence modeling were employed in the 
calculation, and the results of the calculations indicate little difference in 
the prediction of the pressure coefficent. The study of the interaction be- 
tween shock wave and boundary layer in the transonic flow field was not 
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pursued. The mechanism of the transonic shock wave and boundary layer inter- 
actions were reported recently by Roscoe, et al (Ref 38). 

The second flow considered was case 144 in Ref 22. In this case, heat transfer 
effects were included in the calculation. The flow conditions were as follows: 
the inlet Mach number, Ma^ , was 0.16, the inlet total temperature, T^., was 
815°K, the exit Mach number, Ma^, was 0.9, the exit Reynolds number based on 
true chord, Re,,, was 2.43 x 10^, the estimated p -j/ p ex -j t was 1-66. and the free 
stream turbulence intensity, Tu, was 6.5%. The surface temperature distribu- 
tion of the turbine blade is shown in Figure 12. The solid line in the figure 
represents the actual surface temperature used in the calculation and the dot 
symbol represents the data given in Ref 22. The boundary conditions for the 
energy equation, Equation (8), are specified as follows: total temperature is 

held constant at the upstream inlet; surface temperature distribution is speci- 
fied through input by means of Figure 12; conditions on periodic and downstream 
boundaries are treated similarly to other variables at the same boundaries. 

When compared with the Allison experimental data and inviscid predictions in 
Figure 13, the calculations indicate excellent agreement. To indicate the re- 
lationship between the predictions and the experimental scatter, data from 
cases 148 and 158 are also shown in Figure 13. The inviscid pressure calcula- 
tions show close agreement with the computed N-S results. Close agreement was 
expected because this is an on-design case in which viscous displacement ef- 
fects are expected to be small. Significant discrepancies would be expected 
at off-design conditions. Three turbulence models were employed; these 



Figure 12. Surface temperature distribution for case 144. 
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were a mixing length model, a k-l one-equation model, and a k-e two-equa- 
tion model. Computed results based on these three models show little differ- 
ence in the pressure distribution. A vector velocity plot is presented in Fig- 
ure 14. The turning of the flow as it passes through the cascade is evident. 
The flow turning (the exit angle) is predicted by the analysis rather than be- 
ing an input item. The flow acceleration as it passes through the cascade is 
also shown as is the stagnation region and the wall boundary layer development. 

The energy equation was coupled with momentum equations during this calcula- 
tion. In these calculations, the laminar Prandtl number, Pr, was set to 0.73 
and the effective turbulent Prandtl number, Pr^, was set to 0.9. Addition 
of the energy equation to the governing set made little difference to the cal- 
culated pressure distribution, however it allows the calculation of the surface 
heat transfer coefficient. Surface heat transfer prediction on a turbine blade 



Figure 14. Vector plot for C3X cascade, case 144. 
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represents a stringent test due to the flow field character. The turbine pas- 
sage flow field starts from laminar at upstream and then undergoes transition 
and becomes fully turbulent. The calculation of the transitional behavior rep- 
resents a difficult factor in the prediction for the heat transfer coefficient. 
Two cases were considered in employing a mixing length turbulence model; fully 
turbulent flow and transitional flow in the stagnation region. As can be seen 
in Figure 15, the fully turbulent flow case overpredicts the heat transfer co- 
efficient particularly in the stagnation zone. In view of this behavior, a 
transition model was incorporated to investigate the effects on heat transfer. 
Laminar flow was assumed in the region where x/cx is less than 0.2, followed 
by a transitional zone based on the correlation of Dhawan and Narasimha (Ref 
39), and thereafter by fully turbulent flow. The predictions obtained with 
this simple model compare well with the experimental data. 

Although this empirical transition model has given good agreement with data, 
models containing less empiricism would be desirable. One such model, which 
has been used successfully for a variety of transitional flows, is the model 



Figure 15. Comparison of measured and calculated heat transfer 
coefficient distributions of the C3X cascade. 
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of McDonald and Fish (Ref 26), which is based on an integral turbulence energy 
equation and has been used in conjunction with a finite difference boundary 
layer analysis. This model has been generalized from the integral equation 
formulation to the full partial differential equation formulation by Shamroth 
and Gibeling (Ref 28) and incorporated in the N-S code. 

This K-8. model predicts the transport history of the turbulence kinetic en- 
ergy and includes the effect of free-stream turbulence. Although this model 
has not been thoroughly tested under this effort, it may represent a promising 
approach to transitional calculations and is under study by SRA personnel. 

The final 2-D flow considered was case 144 with the film-cooling option. For 
purposes of calculation, it was assumed that air was injected at 30 deg to the 
suction side over 0.8 < x/cx < 0.9 at a velocity of 7% of free stream and the 
local surface temperature of the blade was kept fixed at the same value as the 
nonblowing option. A mixing length turbulence model in conjunction with the 
transition model used in the nonfilm-cooling option for case 144 was employed 
in the calculation. In Figure 16 the distribution of the computed heat trans- 
fer coefficient is depicted for both film-cooling and nonfilm-cooling options. 



Figure 16. The effect of film cooling on heat transfer coefficient 

distributions . 
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The symbol + in the figure indicates the experimental data for the nonfilm- 
cooling case. As can be seen in Figure 16, the heat transfer rate drops to 
nearly zero from the onset of the injection to the trailing edge. This be- 
havior is a consequence of the buffer region of constant temperature cool gas, 
which protects the blade surface from the hotter fluid in the cascade passage. 
The comparison of the pressure distribution for both film-cooling and nonfilm- 
cooling options is shown in Figure 17. The effect of transpiration on the 
pressure distribution is evident. The adverse pressure gradient that is gen- 
erated, the resulting upstream influence, and the subsequent favorable pressure 
gradient that follows it can also be seen in this figure. Temperature, pres- 
sure coefficient, and Mach number contours are given in Figures 18 through 20. 
The effects of the film-cooling are evident. 

CASE 3— 3-D C3X RECTILINEAR TUR8INE CASCADE 

The final case considered was a 3-D demonstration calculation assuming laminar 
flow. The configuration consisted of a C3X cascade situated in the azimuthal- 
radial plane, and bounded in the transverse direction by an endwall and a sym- 
metry plane. The three-dimensionality was introduced by stacking similar 
planes parallel to each other in the direction normal to a fixed endwall (see 
Figure 6). For the calculation, a grid consisting of 100 x 25 x 15 grid points 
in the pseudoazimuthal , pseudoradial , and transverse directions, respectively, 
was constructed. The height of the blade above the endwall (to the symmetry 
plane, midspan) was set to be one axial chord, while the inlet boundary layer 
thickness was 20% of that value. 

The calculation was initiated with an initial condition that consisted of a 
2-D solution with a simple boundary layer correction applied in the vicinity 
of the endwall. The corresponding 2-D solution was obtained with 100 x 25 grid 
points, the inlet Mach number, Ma^ , was 0.15, the inlet Reynolds number based 
on inlet free-stream and time chord, Re 1 , was 730, and the ratio of inlet 
total to downstream exit static pressure, p -t/ p ex -j t » was 1-95. Total temp- 
erature was assumed to be constant. With this initial condition, a converged 
solution was obtained in less than 60 time steps. 
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Figure 17. The effect of film cooling on pressure coefficient 

distributions . 
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b. With film cooling 


TE84-8598 


Figure 19. Pressure contours. 





The computed pressure distributions at different heights above the endwall are 
shown in Figure 21. The pressure side is minimally affected by the endwall, 
remaining at or near the 2-D value run on the same spanwise cross-sectional 
grid, while the suction side, which shows as much as a 15% change over the 2-D 
value near the 30% axial chord location, approaches the 2-D value at 26% span 
above the endwall. Near the endwall, the suction side of the blade is lightly 
loaded compared with the value of the midspan. These differences from the 2-D 
value are due to the effects of secondary flow generated by horseshoe and pas- 
sage vortices. The results are consistent with the expected physics (Ref 6 
and 40 through 43) though different geometry is considered in this effort. 
Figures 22, 23, and 24 show the static pressure contours at endwall, 3.5% span, 
and midspan plane, respectively. The effect of endwall on the pressure distri- 
bution is clear. In Figures 25, 26, and 27, the velocity vector plots are pre- 
sented for the forward portion of the C3X cascade at three different planes 
above the endwall. Near the endwall (within 2.95% spanwise plane) a saddle 
point exists as indicated in the figures. This saddle point moves toward the 
leading edge and disappears beyond the 2.95% spanwise plane. A stagnation 
point forms on the nose of the blade surface beyond the 2.95% spanwise plane. 
These features are consistent with those expected. 
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Figure 21. Three-dimensional rectilinear pressure coefficient 
distribution of the C3X cascade. 
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V. CONCLUSIONS 


The work described in this report has focused on the application of the time- 
dependent ensemble-averaged Navier-Stokes (N-S) equations to transonic turbine 
cascade flow fields. In particular, efforts have focused on an assessment of 
the procedure in conjunction with a suitable turbulence model to calculate 
steady turbine flow fields using an 0-type coordinate system. Three cascade 
configurations have been considered: the two-dimensional (2-0) turbine cas- 

cade of Turner (Ref 21), the 2-D C3X turbine cascade of Hylton et al (Ref 22), 
and the three-dimensional (3-0) C3X rectilinear turbine cascade. The calcula- 
tions were carried out in nonorthogonal body-fitted coordinate systems, while 
grid points were orthogonal near the body surface where viscous flow gradients 
were suitably resolved. A stagnation pressure inflow/static pressure outflow 
boundary condition was employed together with no-slip on the solid surfaces. 

In general, converged solutions were obtained within 80-150 time steps. For 
example, it took 80 time steps to reach a steady-state solution for a 2-D 
laminar Turner turbine cascade flow on a 100 x 30 computational grid starting 
with a uniform velocity (except for a near wall no-slip correction) and a uni- 
form pressure field. The solution was obtained in approximately 570 sec cen- 
tral processing unit (CPU) time of the unvectorized CRAY-1 system. However, 
the code used to implement the algorithm is still a research code with a large 
amount of overhead computation. Recent effort on the code speed-up (Ref 9 and 
44) shows the run time can be reduced by 50% with code restructuring. In addi- 
tion, code vectorization is expected to decrease run time by at least a factor 
of five, leading to a 2-D solution in less than 1 min CPU time. 

Comparisons were made between the predicted and measured surface pressure and 
heat transfer distributions wherever available. In general, the pressure pre- 
dictions were in good agreement with the data. The computed heat transfer re- 
sults also showed good agreement with the data when an empirical transition 
model was used, although further work in the development of laminar-turbulent 
transitional models is indicated. The calculations showed most of the known 
features associated with turbine cascade flow fields. The flow turning, lead- 
ing edge stagnation region, boundary layer development, wake development, film- 
cooling effects, and 3-D boundary layer separation were all clearly observed. 
These results indicate the ability of the present N-S analysis to predict the 
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surface pressure distribution, heat transfer rates, viscous flow development 
for practical turbine cascades operating at realistic flow conditions in 
reasonable amounts of computation time with a suitable turbulence model, and 
plausible boundary conditions. 
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APPENDIX A 
SOLUTION PROCEDURE 


The solution procedure method employs a consistently-split linearized block 
implicit (LBI) algorithm, which has been discussed in detail in Ref 13 
and 34. There are two important elements of this method: 

o the use of a noniterative formal time linearization to produce a fully- 
coupled linear multidimensional scheme that is written in block implicit 
form 

o solution of this linearized coupled scheme using a consistent "splitting" 
[alternating direction implicit (ADI) scheme] patterned after the Douglas- 
Gunn (Ref 35) treatment of scalar ADI schemes. 

The solution procedure method is referred to as a split LBI scheme. The method 
has the following attributes: 

o The noniterative linearization is efficient. 

o The fully-coupled linearized algorithm eliminates instabilities and/or ex- 
tremely slow convergence rates often attributed to methods that employ ad 
hoc decoupling and linearization assumptions to identify nonlinear coeffi- 
cients that are then treated by lag and update techniques. 

o The splitting or ADI technique produces an efficient algorithm that is 
stable for large time steps and also provides a means for convergence ac- 
celeration for further efficiency in computing steady solutions. 

o Intermediate steps of the splitting are consistent with the governing equa- 
tions, and this means that the physical boundary conditions can be used 
for the intermediate solutions. Other splittings that are inconsistent 
can have difficulty in satisfying physical boundary conditions (Ref 34). 

o The convergence rate and overall efficiency of the algorithm are much less 
sensitive to mesh refinement and redistribution than algorithms based on 
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explicit schemes or that employ ad hoc decoupling and linearization assump 
tions. This is important for accuracy and for computing turbulent flows 
with viscous sublayer resolution. 

o The solution procedure method is general and is specifically designed for 
the complex systems of equations that govern multiscale viscous flow in 
complicated geometries. 

The LBI algorithm was later considered by Beam and Warming (Ref 45), but the 
ADI splitting was derived by approximate factorization instead of the Douglas- 
Gunn procedure. Beam and Warming refer to the algorithm as a delta form ap- 
proximate factorization scheme. This scheme replaced an earlier nondelta form 
scheme (Ref 46), which has inconsistent intermediate steps. 

SPATIAL DIFFERENCING AND ARTIFICIAL DISSIPATION 

The spatial differencing procedures are a straightforward adaption of those 
used in Ref 13 and elsewhere. Three-point central difference formulas are used 
for spatial derivatives, including the first-derivative convection and pressure 
gradient terms. This has an advantage over one-sided formulas in flow calcula- 
tions subject to two-point boundary conditions (virtually all viscous or sub- 
sonic flows) in that all boundary conditions enter the algorithm implicitly. 

In practical flow calculations, artificial dissipation is usually needed and 
is added to control high-frequency numerical oscillations that otherwise occur 
with the central-difference formula. 

In the present investigation, artificial (anisotropic) dissipation terms of 
the form 


J J \ax./ ax . 

\ J / J 


(40) 


are added to the right-hand side of each (k-th) component of the momentum equa- 
tion, where for each coordinate direction x^, the artificial diffusivity dj 
is positive and is chosen as the larger of zero and the local quantity 
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y e (° Re flx _1 ) /Re * Here, the local cell Reynolds number Re flx j for the j-th 
direction is defined by 


Re Axj = R® | pu j| (41) 

This treatment lowers the formal accuracy to 0 (Ax), but the functional form 
is such that accuracy in representing physical shear stresses in thin shear 
layers with small normal velocity is not seriously degraded. This latter 
property follows from the anisotropic form of the dissipation and the combina- 
tion of both small normal velocity and small grid spacing in thin shear layers. 

SPLIT LBI ALGORITHM 

Linearization and Time Differencing 

The system of governing equations to be solved consists of either three or four 
equations: continuity and two components of momentum equation in three depend- 

ent variables, p, u, v, and w, or continuity and three components of momentum 
equation in four dependent variables, p, u, v, and w. Using notation similar 
to that in Ref 13, at a single grid point this system of governing equations 
can be written in the following form: 

3H(4>)/8t - D(4>) + S($) (42) 

where <j> is the column-vector of dependent variables, H and S are column-vec- 
tor algebraic functions of 4>, and D is a column vector whose elements are 
the spatial differential operators that generate all spatial derivatives ap- 
pearing in the governing equation associated with each element. 

The solution procedure is based on the following two-level implicit time-dif- 
ference approximations of Equation (42): 


(H n+1 - H n )/At - B(D n+1 + S n+1 ) ( 1-fJ) (D n + S") (43) 
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where, for example, H n+ ^ denotes H(<t> n+ ^) and at = t n+ ^ - t n . The parameter 
B (0.5 < 6 < 1) permits a variable time-centering of the scheme, with a trunca- 
tion error of order [at 2 , (B - 1 / 2 ) at]. 

A local time linearization (Taylor expansion about 4 > n ) of requisite formal 
accuracy is introduced, and this serves to define a linear differential opera- 
tor L (Ref 13) such that 

D n+1 = D° + L n (* n+1 - * n ) + 0( At 2 ) (44) 


Similarly, 


H n+1 = H n + (3H/34>) n (4 > 0+1 - ♦") + 0 (At 2 ) (45) 

s n +i _ s n + ( SS / 3 ^) n ($ n+1 _ 4 , 11 ) + 0 (At 2 ) (46) 

Equations (44 through 46) are inserted into Equation (43) to obtain the follow- 

ing system, which is linear in <t> n+ ^ 


(A - pat L n ) U n+1 - 4 > n ) = at (D n + s n ) (47) 

and is termed an LBI scheme. Here, A denotes a matrix defined by 


A = (3H/H) n ~ Pat (3S/34>) n (48) 

Equation (47) has 0 (at) accuracy unless H = 4 >, in which case the accuracy 
is the same as Equation (43). 

Special Treatment of Diffusive Terms 

The time differencing of diffusive terms is modified to accommodate cross- 
derivative terms and also turbulent viscosity and artificial dissipation coef- 
ficients that depend on the solution variables. Although formal linearization 
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of the convection and pressure gradient terms and the resulting implicit 

coupling of variables is critical to the stability and rapid convergence of 

the algorithm, it does not appear to be important for the turbulent viscosity 

and artificial dissipation coefficients. Because the relationship between 

u and d. and the mean flow variables is not conveniently linearized, 

0 J 

these diffusive coefficients are evaluated explicitly at t during each time 
step. Notational ly, this is equivalent to neglecting terms proportional to 

3y /3<$> or 3d./3$ in l n , which are formally present in Equation 

® J n n 

(44), but retaining all terms proportional to y or d. in both L and 0 . 

° J 

Extensive experience has shown that this has little, if any, affect on the 
performance of the algorithm. This treatment also has the added benefit that 
the turbulence model equations can be decoupled from the system of mean flow 
equations by an appropriate matrix partitioning (Ref 34) and solved separately 
in each step of the ADI solution procedure. This reduces the block size of 
the block tridiagonal systems that must be solved in each step and reduces the 
computational labor. 

In addition, the viscous terms in the present formulation include a number of 
spatial cross-derivative terms. Although it is possible to treat cross-deriva- 
tive terms implicity within the ADI treatment that follows, it is not conveni- 
ent to do so; and consequently, all cross-derivative terms are evaluated ex- 
plicitly at t n . For a scalar model equation representing combined convection 
and diffusion, it has been shown by Beam and Warming (Ref 47) that the expli- 
city treatment of cross-derivative terms does not degrade the unconditional 
stability of the present algorithm. To preserve notational simplicity, it is 
understood that all cross-derivative terms appearing in L n are neglected but 
are retained in D n . Neglecting terms in L n has no effect on steady solu- 
tions of Equation (47), because 4>° +1 - 4>° = 0. and thus Equation (47) 
reduces to the steady form of the equations: D n + S n = 0. Aside from 

stability considerations, the only effort of neglecting terms in L n is to 
introduce an 0 (At) truncation error. 
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Consistent Splitting of the LBI Scheme 


To obtain an efficient algorithm, the linearized system Equation (47) is split 
using ADI techniques. To obtain the split scheme, the multidimensional opera- 
tor L is rewritten as the sum of three one-dimensional (1-D) suboperators 
(i = 1, 2, 3) each of which contains all terms having derivatives with respect 
to the i-th coordinate. The split form of Equation (47) can be derived either 
as in Ref 13 and 34 by following the procedure described by Douglas and Gunn 
(Ref 35) in their generalization and unification of scalar ADI schemes, or us- 
ing approximate factorization. For this system of equations, the split algo- 
rithm is given by 


(A - BAtL*) ($* - *") = At (D n + S n ) (49a) 

(A - BAtL”) (♦** - - A (♦* - ♦") (49b) 

(A - BAtL") ($ n+1 - *“) - A (♦** - ♦”) (49c) 

where <t>* and $** are consistent intermediate solutions. If spatial deriva- 
tives appearing in L.. and D are replaced by three-point difference formulas, 
as indicated previously, then each step in Equation (49) can be solved by a 
block-tridiagonal elimination. 

Combining Equations (49a-c) gives 

(A - BAtL°) A _1 (A - BAtL") A _1 (A - BAtL") (♦ n+1 - ♦") - At (D n + S n ) (50) 

2 

which approximates the unsplit scheme Equation (47) to 0 (At ). Because 
the intermediate steps are also consistent approximations for Equation (47), 
physical boundary conditions can be used for and <|>** (Ref 13 and 34). 
Finally, since the suboperators L.. are homogeneous operators, it follows from 
Equations (49a-c) that steady solutions have the property of <|> n+1 = $* = <$>** 

= <$> n and satisfy 


n 


D + S 


n 


(51) 
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The steady solution depends only on the spatial difference approximations used 
for Equation (51) and does not depend on the solution algorithm itself. 
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APPENDIX B 
USER'S MANUAL 

The user's manual that follows describes two separate computer codes that are 
used to obtain solutions for cascade flow problems. A nonorthogonal coordinate 
system is generated through COORD so that the physical boundaries and periodic 
boundaries coincide with coordinate surfaces. The output from COORD consists 
of results printed in the output file and results written on a temporary file. 
This latter set of results must be saved and used as input for the second pro- 
gram, MINT. 

Figure 28 shows the overall program flow of COORD. Table I describes the COORD 
namelist input. The MINT code combines a block data program (BLKDAT) contain- 
ing pertinent data statements, a main program (DAL) and a series of subroutines 
to perform the required calculations. Figure 29 shows the overall program 
flow. Figure 30 illustrates the input and initialization procedures. Table II 
describes the MINT namelist input, and Table III lists the major FORTRAN vari- 
ables in MINT. Figure 31 provides a global description of execution control. 
These program flow charts only provide a means of guidance for interested users 
wishing to consult the program listing. 

SAMPLE INPUT CARDS AND PRINTED OUTPUT 

A sample two-dimensional (2-D) case was run to illustrate the set up of input 

paramaters and typical printouts of the COORD program and MINT code, respec- 

tively. Table IV lists sample inputs for COORD. For the program, all param- 
eters in the namelist input &INPT and 8.DATA1 must be specified. Two sequential 
&DATA1 are required in the namelist input. The first set is used for the 
parameterization of the inner loop and the second set is used for that of the 
outer loop. Sample input for the MINT code is listed in Table V. For the MINT 
code, the namelist data shown is for a new calculation. In the case of a re- 
started calculation the input is identical except for the restart flag IREST 

set in &READ1, and the variables NT and DT in &READ9. It is not advisable to 

change any other values on a restart. 

In COORD, the code output first prints out the input airfoil coordinate data, 
as shown in Table VI. These geometry data are followed by the print-out of 
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the namelist INPT and two sequential DATA1 . Finally, the grid coordinates' 
associated metric data, coordinate intersection angles, and Jacobian are 
printed out. The grid coordinates are written by means of a binary write on 
TAPE9 and must be saved for input to the MINT code. 

In the MINT code, the code output first prints out a series of dimensionless 
parameters (DIM1-DIM10, DIM12, and DIM14) and the dimensionless, referenced 
total temperature, total pressure, and total enthalpy as shown in Table VII. 
This output is followed by the finite difference coefficients for first and 
second derivatives in both directions. In each direction, three lines are 
written. The first and third lines give one-sided difference weights at the 
lower and upper boundaries; the second line of each set gives central differ- 
ences used for the interior points. Six numbers are written on each line; the 
first set of three values represents the first derivate coefficients and the 
second set of three values represents the second derivative coefficients. 

The next output item is the namelist data &READ1 and &READ9. The values of 
some parameters may be changed by the internal operations. The following items 
are the print-outs of the local surface static pressure, boundary layer thick- 
ness, and heat transfer coefficients around the blade surface, respectively. 

The upstream total pressure is also indicated in the print-out. The value of 
total pressure, PTOTI, will be reached in IDTAB (through namelist input) time 
steps if the PTOTI is greater than the reference total pressure. A summary 
print, which gives not only the maximum relative change over a time step by 
SSTES1 but also the location, is written at each time step. The maximum rela- 
tive changes of each dependent variable are also given. RESMAX is the maximum 
residual of the equations solved and indicates how well the steady-state equa- 
tion is being satisfied. 

Following the summary data of the final time step of the run, the pressure co- 
efficient and the heat transfer coefficient are written against the blade sur- 
face coordinates. The distribution of the heat transfer coefficient along the 
blade surface is zero for the case of IHSTAG=2 (i.e. constant stagnation 
enthalpy option). Finally, a complete flow field print-out is produced. 
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Figure 28. The overall program flow for program COORD 










Figure 29. The overall program flow for program DAL 





























Namelist or 
variable name 

&INPT 

NRAD 

NTHETA 

RMAXO 

RMAX1 

RMAX2 

DAMP 

NE 

LED 

DMX(I) 

DMY(I) 

XM2, XM3 
YM2, YM3 

EXTF , EXTR 

&DATA1 

NCLUST 

CLPY(I), I = 1 
CLPX(I) , I = 1 
ETAP(J), I = 1 
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Table I. 

COORD namelist input description. 

Description 

number of pseudoradial points 
number of pseudoazimuthal points 
cascade passage spacing/2 

normal distance between blade surface and inner-second- 
ary loop 

normal distance between outer loop and outer-secondary 
loop 

damping parameter in pseudoradial direction 
number of airfoil input data points 
airfoil data point number at airfoil leading edge 
array for x-location of airfoil input data points 
array for y-location of airfoil input data points 
input outer loop location (see Figure 32) 

extension of periodic line (see Figure 32) 


total number of the interior cluster points [A cluster 
point is the sequential number of the selected grid 
point that must coincide with a particular predetermined 
value of the z-coordinate. Accordingly, pairs of (1, 
YFIRST) and (NUMDZ+2, YLAST) are also cluster points, 
but they are boundary cluster points.] 

- 60 the z-coordinate of the cluster point (Both boundary 

and internal cluster points must be specified.) 

- 60 the sequential number of the grid point corresponding 

to CLPY(I) 

- 60 the sequential number of the grid point defined as the 

pivot point (The grid spacing will have the fastest 
variation at a pivot point. For each of the interior 
cluster points there will be a pair of pivot points: 
one ahead of the cluster point, the other after the 
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Table I. (cont) 

Namelist or 

variable name Description 

cluster point. However, only one pivot point will be 
associated with each of the boundary cluster points.) 

ALPH(J), 1=1 - 60 width parameter specifying width (in terms of the number 

of grid points) in which 90% of the grid size variation 
takes place around the pivot point ETAP(J) 

> 0 decreasing grid size 

< 0 increasing grid size 

NEND =0 no stretching at YFIRST and YLAST 

= 1 stretching at YFIRST only 

= 2 stretching at YLAST only 

= 3 stretching at YFIRST and YLAST 

RHO approximate ratio of the grid size at CLPY(2) to the 

(always > 0.0) maximum grid size in the interval CLPY(l) < Z < 

CLPY( 2) 

no stretching at YFIRST (used with NEND = 0 or 2) 

40 approximate ratio of the grid size at CLPY(K+1) to the 
maximum grid size in the interval CLPY(K) < Z < 

CLPY(K+1 ) 

=1.0 no grid variation at CLPY(K+1) 

BETA0 calculated (It indicates the first derivative of the 

(always > 0.0) z-coordinate with respect to the computational 

coordinate at YFIRST.) 

BETA(L), L = 1 - 60 calculated [It indicates the ratio between the grid 
(always > 0.0) sizes on both sides of the pivot point ETAP(L+1).] 


Notes: 

The airfoil x-y data points DMX and DMY must be input in order around the air- 
foil starting from any arbitrary location on the airfoil. The computer code 
then nondimensionalizes the coordinates with respect to the axial chord and 
places the origin at the location of maximum x (near the trailing edge) with 
the airfoil oriented as shown in Figure 32. The units of DMX and DMY are 
arbitrary while RMAX0, RMAX1 and RMAX2 are nondimensional . 

After this is done, the grid point numbering is reordered so that the first 
point terminates in the trailing cap, e.g., line AA' in Figure 2. All pseudo- 
azimuthal grid point locations are referenced to this line AA'. Hence, CLPX 
and CLPY refer to the grid point locations and the azimuthal parameterized co- 
ordinate value, respectively. 


= 1.0 

RATIO(K) , K = 1 - 
(always > 0.0) 


71 



y 



Figure 32. Cascade geometry. 


Namelist or 
variable name 


Table II. 

MINT namelist input description. 


Description 


&READ1 

IREST = 0 
= 1 

IOTAPE = 10 
INTAPE = 9 


A new calculation is being started. 

The case is being run from the restart files. 

output unit number for the dependent variable array 
restart data 

input unit number for the dependent variable array 
restart data 


I0TAP1 = 20 
INTAP1 = 19 
&READ9 
NUMDX 

NUMDY 

NUMDZ 

TWOD = T 
= F 

CLENG 

WREF 

DENSR 

TREF 

VISCR 

AVISC(IDIR.IEQ) 


output unit number for the namelist restart data 
input unit number for the namelist restart data 


number of interior grid points in the pseudoazimuthal 
direction (x or y 1 direction); total number of points 
in this direction = NUMDX + 2, NUMDX < 128 

number of interior grid points in the spanwise 
direction (3D only) 

number of interior grid points in the pseudoradial 
direction (z or y 3 direction); total number of points 
in this direction = NUMDZ + 2, NUMDZ < 28 

logical variable for 2-D calculation 

logical variable for 3-D calculation 


reference length — m 

reference velocity — m/sec 

reference density — kg/m 3 

reference temperature — °K 

reference dynamic viscosity — kg/m sec 

artificial dissipation parameter a 

IEQ = 1-5 and IDIR = 1-3 

IEQ = 1 indicates the x-momentum equation. 

IEQ = 2 indicates the y-momentum equation. 

IEQ = 3 indicates the z-momentum equation. 

IEQ = 4 indicates the continuity equation. 
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Namelist or 
variable name 


Table II. (cont) 


NT 

OT 


DTMIN 

DTMAX 

ITEST 

IPRINT 

IVARPR(IV) 

= 0 
= 1 


Description 

IEQ = 5 indicates the energy equation. 

IEQ = 16 indicates the turbulent kinetic energy. 

IEQ = 17 indicates the turbulent dissipation energy rate. 
IDIR = 1 indicates the x second derivative term. 

IDIR = 2 indicates the y second derivative term. 

IDIR = 3 indicates the z second derivative term. 

For example, AVISC(3, 1) is the value of o used for 

the artificial dissipation term a^U/aZ^ in the 
x-momentum equation. Even if the y-momentum equation is 
not solved, the corresponding AVISC must be supplied. 
Default values are 0.0; recommended values are 0.50 
initially, followed by runs at 0.2. 

number of time steps to be run 

initial nondimensional time step (If DT is omitted on a 
restart, DT will be set to the value it was at the 
termination of the last run.) 

minimum nondimensional time step for this run 

maximum nondimensional time step for this run 

Steady-state test is performed every ITEST time steps. 
(Default value is 1 .) 

Complete flow field print-outs are provided every IPRINT 
time steps. 

optional print control flag for variable IV 
suppress print-out of variable IV 
normal print-out of variable IV 

IV IVARPR(IV) 

1 transverse velocity, u 

2 spanwise velocity, v 

3 streamwise velocity, w 

4 density, p 

5 enthalpy, h 

16 turbulent kinetic energy, k 

17 turbulent dissipation rate, e 

26 pressure, P 

27 temperature, T 

28 effective viscosity, vef f 

33 mixing length, 9. 
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Table II. (cont) 


Namelist or 
variable name 

IDUMP1 

= 1 
= 2 

I PLOT 

= 0 
> 0 

JDBLU 

ODBLl 

KDBLU 

KDBLL 

IHSTAG 

= 0 
= 2 

FLWNG 

TGAS 


Description 

35 Mach number, Ma 

36 total pressure, P-j- 

print initial flow field for this run 
no initial print-out 


no plot file (TAPE1) written 

plot file written at time step increment IPLOT 

point number at the end of the periodic line in the 
leading edge of the cascade (see Figure 32) referenced 
to the grid point corresponding to line AA 1 in Figure 2 

point number at the end of the periodic line in the 
leading edge of the cascade (see Figure 32) referenced 
to the grid point corresponding to line AA' in Figure 2 

point number at the end of the periodic line in the 
trailing edge of the cascade (see Figure 32) referenced 
to the grid point corresponding to line AA 1 in Figure 2 

point number at the end of the periodic line in the 
trailing edge of the cascade (see Figure 32) referenced 
to the grid point corresponding to line AA 1 in Figure 2 


static enthalpy option 

constant stagnation enthalpy option 

exit flow angle used for defining the wake region and 
boundary extrapolation, referenced to streamwise axis, 
x, as shown in Figure 32 

total temperature of the incoming gas at the inlet 
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Namelist or 
variable name 


Table II. (cont) 


TWALL1 

IBLOW 

= 0 
= 1 

BLOU 

BLOW 

PTOTI 

IVISC 

= 1 
= 2 


Description 

surface temperature distribution along the blade 
surface (nondimensional ) , referenced to pseudoazimuthal 
grid point location around body 


without film-cooling option 
with film-cooling option 

U velocity component of the blowing velocity profile 
for the film-cooling option (nondimensional), 
referenced to pseudoazimuthal grid point location 
around body 

W velocity component of the blowing velocity profile 
for the film-cooling option (nondimensional), 
referenced to pseudoazimuthal grid point location 
around body 

total pressure at the inlet of the cascade flow field 
(nondimensional ) 

index for turbulent modeling 

laminar flow constant viscosity 

laminar flow Sutherland's law 


= 3 

= 4 
= 5 


IOTA 

IDTAB 


mixing length model 
k-a, one-equation model 
k-c two-equation model 

the time step number at which the inlet total pressure 
begins to be increased linearly (default IDTA = 0) 

The number of time steps used to increase inlet total 
pressure up to PTOTI begins at IDTA (default IDTAB = 
10 ). 


Notes: 

The nondimensional time is defined as CLENG/WREF. The recommended values of 
the time steps can be obtained from the sample run in this User's Manual. 
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Table III. 

List of major FORTRAN variables in HINT. 


FORTRAN 

symbol 

Common 

block 

Description 

AC( I , J , K) 

BLK1 

dependent variable array 

ACG( J, J) 

BLK1 

geometry data array 

AN( I , J) 

BLKM 

array storing time term linearized coefficients 

APR( I , J) 

PRNT 

print output array 

AVANDR 

TURB 

damping constant 

AVISC(I.J) 

MISC2 

artificial dissipation parameter 

C( I , J , K) 

BLKM 

coupled matrix array storage 

CLENG 

CREF 

reference length 

CMACH 

MISC2 

reference Mach number 

0 

VARNO 

index for divergence 

01(1. J.K) 

BLKM 

array storing first sweep linearized 
coefficients 

D2( I , J ,K) 

BLKM 

array storing second sweep linearized 
coefficients 

D3( I , J.K) 

BLKM 

array storing third sweep linearized 
coefficients 

DENSR 

CREF 

reference density 

DFW( I , J ,K) 

ADI 7 

difference weight array 

DIM1 

NONO 

inverse Reynolds number 

DIM2 

NONO 

reference pressure/reference dynamic head 

0IM3 

NONO 

reference pressure/( reference density x 
reference enthalpy) 

DIM4 

NOND 

1 .0/(Rey x P r ) 

DIM12 

NOND 

2.0 x 0IM1 

OS 

VARNO 

index for dissipation 

DT 

MISC2 

time step 
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Table III. (cont) 


FORTRAN 

svmbol 

Common 

block 

Description 

OTCON 

MISC2 

inverse step 

DTMAX 

M I SC 2 

maximum allowable time step 

OTMIN 

MISC2 

minimum allowable time step 

E(I,J.K) 

BLKM 

coupled matrix array storage 

GRIO(I) 

GTRAN 

grid distribution parameter 

H 

VARNO 

index for enthalpy 

11 

MGAUS 

lower limit for matrix inversion 

IADI 

AD 1 1 

ADI sweep number 

IBC 

AD 1 1 

boundary condition boundary parameter 

IDT 

MISC2 

time step index 

IDTADJ 

MISC2 

time step control parameter 

IDUPM1 

OUTA 

parameter controlling initial station print 
(see &READ9) 

IEQ 

AD 1 1 

equation number 

IGPRT(I) 

GEOl 

geometry print control 

IL 

MGAUS 

upper limit for matrix inversion 

IPRINT 

MISC2 

print interval parameter (see &READ9) 

IRERUN 

MISC2 

restart write control parameter 

IREST 

MISC2 

restart read control parameter (see &READ1) 

IVARPR(I) 

MISC2 

print parameter (see &READ9) 

JADI 

AD 1 1 

ADI sweep parameter 

JEQBC( I , J,K) 

ADI1 

boundary condition type parameter 

JX 

ADI2 

direction 1 grid point index 

KZ 

ADI2 

direction 3 grid point index 

LX 

ADI2 

direction 1 grid point index 
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Table III. (cont) 


FORTRAN 

symbol 

Common 

block 

Description 

LX1 

ADI3 

first direction 1 interior point 

LX2 

ADI3 

last direction 1 interior point 

LY 

A0I2 

direction 2 grid point index 

LY1 

ADI3 

first direction 2 interior point 

LY2 

ADI3 

last direction 2 interior point 

LZ 

ADI2 

direction 3 grid point index 

LZ1 

ADI3 

first direction 3 interior point 

LZ2 

ADI3 

last direction 3 interior point 

LZPRNT 

MISC2 

3-D print control 

MEQS 

ADI1 

number of equations to be solved 

NT 

MISC2 

number of time steps to be run 

NUMDX 

MISC2 

number of interior direction 1 points 

NUMDY 

MISC2 

number of interior direction 2 points 

NUMDZ 

MISC2 

number of interior direction 3 points 

NX1 

ADI4 

first grid point — direction 1 

NX2 

ADI4 

last grid point — direction 1 

NY1 

A0I4 

first grid point— direction 2 

NY2 

ADI4 

last grid point — direction 2 

NZ1 

ADI4 

first grid point — direction 3 

NZ2 

ADI4 

last grid point — direction 3 

P 

VARNO 

index for pressure 

PCNT1 

MISC2 

time step control parameter 

PCNT2 

MISC2 

time step control parameter 

PREF 

CREF 

reference pressure 

PRNDL 

CREF 

Prandtl number 
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Table III. (cont) 


FORTRAN 

symbol 

PTOT 

R 

REY 

SN(I) 

SSTEST 

T 

TAUW 

TREF 

TTIME 

TTOT 

U 

USTAR 

V 

VISCL 

VISCR 

VS 

W 

WREF 

XGMAX(I) 

XGMIN(I) 

XPLUS 


Common 


block 


Description 

BCCON 

total 

pressure 

VARNO 

index 

for density 

CREF 

Reynolds number 


BLKM 

array storing source term linearized 

coefficient 

MISC2 

maximum change in variable across the 

time step 

VARNO 

index for temperature 


TURB 

wall shear 


CREF 

reference temperature 


MISC2 

cumulative time 


BCCON 

total temperature 


VARNO 

index for direction 1 velocity 


TURB 

dimensionless velocity 


VARNO 

index for direction 2 velocity 


TURB 

laminar reference viscosity 


CREF 

reference viscosity 


VARNO 

index for viscosity 


VARNO 

index for direction 3 velocity 


CREF 

reference velocity 


GRID1 

maximum coordinate value 


GRID1 

minimum coordinate value 


TURB 

dimensionless distance from surface 
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Table VII. (cont) 


oo ssssssg^sssssss£2sf|«||s|ss||«||os||ss|o|ov? ! 

I ^r^rr!r^rvmr^ 

0’*>tO^* > 'T i, r ,>r> ' r " r t . 'oopooo • OOOOOOOOOOOOOOOOOOO 

^oo«oo«ooooooooooooooooooooooooooooooooo o oooooooooocoo^gooooooooooo 

ooooooooocoooooooogcoeooooooooogoo-oo-o^ c ° ?t?“ ~ • *°* • • • • : 

? oooooooac,ooo # *♦...♦ • • . • * -I*:**:: s ?s^^s--ssiSj;^s3Sss;s;5KS 

sssksss:™;::^^^^ s ^ww522222222222^^- 

• * ’ ’ 'oooooOOOOOOOO . oooooooooooooooc oo 

! ISgPilliSI 

““S^SS^^S£SSjg?^j^H5HHs2I^2SSS-;St:S||oS?SS^; -1 ;SSS!^SSS:£SSS2SrifS5fS2?SS£ 
SESS£:2SSSSS^rr =5 ; S -«^-22SS5S52S2Sr2SS- s ~o S; -- ° 

ssssssss»5«5-ssss^ w 

^OOOOOOOO^OOOOOOOOOOOOOOOOOOOOOOOOOO^ 00000000=0000000000000=0000 

sssssssssssssssfS2Sssssssggg5gf?22^^Ts?T?22? • ^^SSSiaasas^ 

S;sHs^5s^^sK^s53SSS^“ ?sss:: 5i^S2 'I Ss£^^s~K5s^.™«* 

°. 000000000=00000000000000000 

ooooooooooooooooooooooooooooooooooooo^ooooooooooo s 000 000000000000000000000000 

0000000000000000000000000000000=00-00^^ ^^ O oooooooooo _ ........ ♦ ♦ * *ii*,*, 

,1 sS^Ss^^sJi 

° r*rgrwc 4 r*f*cvr\tf* rvr\*r.<rg r*r«r<Kii ''.*'^«0 ^ * * 

sjj^)»«*« ( '"* oira; ' i “ o * * * -oooQtjoao 

* «~«aoooooooo • ooo o oooo oooo oooc. o woo 

0000=00000000=000000000000000000000000 S ooooooooooooooooooooooooooo 

ooooooooooooooo gsssssssg °sgsgsgssssoSoSoooooooooo o ^-“^-^nrrrr?. ...♦♦•• 

ii^ii^^^“““^^s^“ssssssgs5s :; =-- s:s s s ss!:ssns --*3 = Essssi££^£SS'.sc^Sii^' 

KHssr.s^^5S?.s^s^s^r.» 

° ^00000000=0000000000000000 

^oooooooooooooooooooooooooooooooooo^oooo „ 0 ooooooooooooo SgS2gSSS ° SSSS 

•SSS2Soooo|ooooooooooooooooooooooo^ gg 2?T????T?T2TT £ 

1 SHEkskksssss^ 

- Ooooooooooooooooooooooooooo 

1 SSSSKSSKKWS8S 

• SoSoSSoSoSSoooooooooooooo«=»oo«=c»o^oooooo^^^^^^^^^ ♦ * 

2S En?fes?wss=«K«^s^ 

s -T^"n r /:^ r ^22ri^iri^ooooo 

^ooooooooooooooooooooooooooooo | SSSSSSSSSSSSSSSSSSSSSSSSffS 

ssoooooooooSoooooooo i 

“SS^SKK^S^SSS^^ ~f W&^»«SSS=SS5KKS»2 

^S^kksSss25^SHSSSSsS32=s?wss i ^.rr.^r:^r:rr.-.r:-rr.';-'. 

• ’ ' ‘^oeooooooo ° oooooooooooo oooo ooo oooooooc 

! !i«llll 

rss5sssss2ssussssgsss|||s|ss||||0||ss|^ EssssSssSSSsssscasss™ 

s 222222 25 22222 ZZ&i 

• • • O^oooo OOOOOOOOOOOO O OOO OOO o 

OOODOCJOOOOOOOOOOOOtrvJt’--' 


108 



Table VII. (cont) 
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Table VII. (cont) 
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Table VII. (cont) 
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APPENDIX C 

LIST OF ACRONYMS, ABBREVIATIONS, AND SYMBOLS 


Term 

Definition 


A 4 

van Driest damping coefficient 


ADI 

alternating direction implicit 


ADVR 

axial velocity density ratio 


b 

total or effective viscosity 


c p 

specific heat at constant pressure 


D 

damping parameter 


d 

distance to the nearest wall 


d 4 

dimensionless distance to the nearest wall 


D 

deformation tensor 


h 

enthalpy 


h/h 0 

normalized heat transfer coefficient 


hi, h 2 , h 3 

normal distance between loops 


I 

identity tensor 


k 

turbulence kinetic energy 


k B 

bulk viscosity coefficient 


L8I 

linearized block implicit 


a 

mixing length 


loo 

mixing length in the core flow region 


Ma 

Mach number 


MINT code 

multidimensional, implicit, nonlinear, time-dependent 
code 


N-S 

Navier-Stokes 


P 

pressure 


Pr 

1 

laminar Prandtl number 


P(r,s) 

position vector 
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Term 


Prt 

Ps 

P S /P T 

Pt 

Q 

q 

Qt 

R 

r 

Re 

Rt 

s 

SRA 

T 

t 

To 

Tt 

T w /Tg 

u 

U* 

u 

— * 
u 

u e 


Definition 

effective turbulent Prandtl number 
static pressure 

static pressure/total pressure 
total pressure 
mean heat flux vector 
magnitude of the velocity 
turbulent heat flux vector 
universal gas constant 
radial parameter 
Reynolds number 

local ratio of turbulent to laminar viscosity 

arc length parameter 

Scientific Research Associates 

temperature 

time 

reference temperature 
stagnation temperature 
wall to gas temperature ratio 
velocity 
velocity vector 

velocity component in x-direction 
fluctuation velocity vector 
edge velocity 
friction velocity 
velocity component in y-direction 
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Term 


Definition 


w 

velocity component in z-direction 

w e 

w at the edge of the boundary layer 

X, X] 

Cartesian coordinate in the transverse direction 

x/cx 

fraction of axial chord 

x te 

trailing edge location 

y, x 2 

Cartesian coordinate in the spanwise direction 

y 1 , y 2 , y 3 

computational coordinates 

z. * 3 

Cartesian coordinate in the streamwise direction 

l-D 

one-dimensional 

2-D 

two-dimensional 

3-D 

three-dimensional 

AX 

grid spacing 

<5 

boundary layer thickness 

4 PS 

pressure surface trailing edge boundary layer thickness 

4 SS 

suction surface trailing edge boundary layer thickness 

c 

turbulence energy dissipation rate 

K 

von Karman constant 

K 

mean thermal conductivity 


dynamic viscosity 


turbulent viscosity 

v a rt 

artificial dissipation 

n, f 

computational coordinates 

if 

molecular stress tensor 

*T 

turbulent stress tensor 

P 

density 
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Term 


Definition 


0 

T 

n 

T xx* 

* 


artificial dissipation parameter 
time 

local shear stress 

x etc. component of stress tensor 

xy * 

mean flow dissipation rate 
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